telluric: manipulación interactiva de datos geoespaciales en Jupyter con Python

Juan Luis Cano - 2018-06-07 @ Jornadas de SIG Libre

https://github.com/satellogic/telluric-talks

telluric

1. Introducción

Sobre Satellogic 🌍

  • Estamos construyendo una constelación de satélites de observación de la Tierra
  • Estamos creando infraestructura para recopilar y almacenar datos (fundamentalmente imágenes) de esta constelación
  • Estamos extrayendo información de estos datos así como de fuentes externas
  • ¡Estamos contratando!

https://www.satellogic.com/

Satellogic

Sobre mí

  • Ingeniero Aeronáutico especializado en mecánica orbital 🛰
  • Socio fundador y presidente de la asociación sin ánimo de lucro Python España y organizador de la PyConES 🐍
  • Software Engineer en el equipo encargado de recopilar, almacenar y acceder a datos geoespaciales en Satellogic 🌍
  • Defensor del Software Libre y entusiasta de Python para ciencia e ingeniería

Juan Luis Cano

2. Software existente

  • Python es un lenguaje de programación de propósito general con presencia en el mundo GIS
  • Además de lenguaje de scripting (QGIS, ArcGIS) existe un ecosistema de bibliotecas para uso geoespacial
  • Acceso a bases de datos geoespaciales, entrada/salida, procesamiento, visualización...
  • Sin embargo, este ecosistema es un poco complicado de descifrar:

Landscape

Resumen de proyectos clave:

Python C/C++
pyproj PROJ.4
Fiona OGR
Shapely GEOS
rasterio GDAL

Ventajas:

  • API "pythonica"
  • "Filosofía UNIX" (hacer una sola cosa y hacerla bien)
  • Bibliotecas con un alcance específico, ligeras y reutilizables

Desventajas:

  • La integración corre a cargo del usuario
  • El manejo de proyecciones, en algunos casos, también

3. telluric

telluric is an open source (MIT) Python library to manage vector and raster geospatial data in an interactive and easy way.

Importar para uso interactivo:

In [1]:
import telluric as tl
from telluric.constants import WGS84_CRS, WEB_MERCATOR_CRS

Basic geometry definition using GeoVector: in the simplest case, the bounds and the projection (CRS)

In [2]:
ebro = tl.GeoVector.from_bounds(
    xmin=0, ymin=40, xmax=1, ymax=41, crs=WGS84_CRS
)
print(ebro)
GeoVector(shape=POLYGON ((0 40, 0 41, 1 41, 1 40, 0 40)), crs=CRS({'init': 'epsg:4326'}))
In [3]:
ebro
Out[3]:

También se pueden usar geometrías de Shapely:

In [4]:
from shapely.geometry import Polygon

ebro2 = tl.GeoVector(
    Polygon([(0, 40), (1, 40.1), (1, 41), (-0.5, 40.5), (0, 40)]),
    WGS84_CRS
)
ebro2
Out[4]:

Propiedades geométricas

In [5]:
ebro.centroid
Out[5]:
In [6]:
ebro.area  # Real area in square meters
Out[6]:
9422706289.175217

Relaciones geométricas

In [7]:
ebro.within(ebro2)
Out[7]:
False
In [8]:
print(ebro.difference(ebro2))
GeoVector(shape=MULTIPOLYGON (((0 40.66666666666666, 0 41, 1 41, 0 40.66666666666666)), ((1 40.1, 1 40, 0 40, 1 40.1))), crs=CRS({'init': 'epsg:4326'}))
In [9]:
ebro.difference(ebro2)
Out[9]:
  • Features: GeoVector + atributos
  • FeatureCollections: una secuencia de features
In [10]:
gf1 = tl.GeoFeature(
    ebro,
    {'name': 'Delta del Ebro'}
)
gf2 = tl.GeoFeature(
    ebro2,
    {'name': 'Delta del Ebro + Castellón'}
)
print(gf1)
print(gf2)
GeoFeature(Polygon, {'name': 'Delta del Ebro'})
GeoFeature(Polygon, {'name': 'Delta del Ebro + Castellón'})
In [11]:
fc = tl.FeatureCollection([gf1, gf2])
fc
Out[11]:

Datos raster

In [12]:
# This will only save the URL in memory
rs = tl.GeoRaster2.open(
    "https://github.com/mapbox/rasterio/raw/master/tests/data/rgb_deflate.tif"
)

# These calls will fetch some GeoTIFF metadata
# without reading the whole image
print(rs.crs)
print(rs.band_names)
CRS({'init': 'epsg:32618'})
[0, 1, 2]
In [13]:
rs.footprint()  # Bounding box
Out[13]:
In [14]:
rs
Out[14]:
In [15]:
rs.image[:,200:300, 200:240]
Out[15]:
masked_array(
  data=[[[15, 20, 23, ..., 41, 38, 29],
         [17, 15, 23, ..., 41, 46, 46],
         [18, 17, 17, ..., 39, 39, 44],
         ...,
         [106, 49, 62, ..., 46, 58, 50],
         [255, 255, 255, ..., 38, 37, 43],
         [255, 229, 255, ..., 40, 31, 35]],

        [[94, 98, 100, ..., 129, 98, 62],
         [98, 96, 100, ..., 127, 127, 133],
         [98, 96, 94, ..., 126, 126, 129],
         ...,
         [145, 79, 116, ..., 47, 58, 51],
         [255, 255, 255, ..., 43, 41, 45],
         [255, 249, 255, ..., 41, 33, 37]],

        [[131, 131, 131, ..., 161, 125, 86],
         [133, 133, 135, ..., 164, 160, 161],
         [133, 132, 129, ..., 159, 160, 157],
         ...,
         [141, 77, 108, ..., 29, 34, 31],
         [255, 255, 255, ..., 31, 30, 31],
         [255, 255, 255, ..., 23, 19, 27]]],
  mask=[[[False, False, False, ..., False, False, False],
         [False, False, False, ..., False, False, False],
         [False, False, False, ..., False, False, False],
         ...,
         [False, False, False, ..., False, False, False],
         [False, False, False, ..., False, False, False],
         [False, False, False, ..., False, False, False]],

        [[False, False, False, ..., False, False, False],
         [False, False, False, ..., False, False, False],
         [False, False, False, ..., False, False, False],
         ...,
         [False, False, False, ..., False, False, False],
         [False, False, False, ..., False, False, False],
         [False, False, False, ..., False, False, False]],

        [[False, False, False, ..., False, False, False],
         [False, False, False, ..., False, False, False],
         [False, False, False, ..., False, False, False],
         ...,
         [False, False, False, ..., False, False, False],
         [False, False, False, ..., False, False, False],
         [False, False, False, ..., False, False, False]]],
  fill_value=999999,
  dtype=uint8)
In [16]:
# Crop fetch only the data required to build the raster
rs.crop(rs.footprint().buffer(-50000))
Out[16]:
In [17]:
rs[200:300, 200:240] # the image is streched due to the presentation layer
Out[17]:
In [18]:
# Rasterizing Vectors
# Geting OpenStreetMap data of Tel Aviv
fc = tl.FileCollection.open("data/telaviv.geojson")
fc
Out[18]:
In [19]:
streets_raster = fc.rasterize(dest_resolution=1)
streets_raster
Out[19]:

Computación distribuida

  • Podemos usar .crop y .rasterization para distribuir la carga en varias máquinas
  • El formato Cloud Optimized GeoTIFF (COG) permite optimizar las lecturas parciales de datos raster
  • Podemos almacenar los datos en un blob storage (Azure, AWS S3)
In [20]:
from telluric.vectors import generate_tile_coordinates

list_of_regions = list(generate_tile_coordinates(rs.footprint(), num_tiles=(10,10)))

# This line is to visualize and not required  
tl.FeatureCollection.from_geovectors(list_of_regions)
Out[20]:
In [21]:
chunk0 = rs.crop(list_of_regions[17])
chunk0
Out[21]:
In [22]:
import dask.multiprocessing
from dask.delayed import delayed
from dask.base import compute

def worker(raster_filename, roi):
    raster = tl.GeoRaster2.open(raster_filename) #Lazy loading of the raster
    chunk0 = raster.crop(roi)
    return chunk0.image.max()

RASTER_URL="https://github.com/mapbox/rasterio/raw/master/tests/data/rgb_deflate.tif"
items = [delayed(worker)(RASTER_URL, roi) for roi in list_of_regions]

maxs = compute(*items, get=dask.multiprocessing.get)

4. Uso interno en Satellogic

  • Visualizar capturas y trazas de nuestros satélites
  • Estimar la cobertura y hacer cálculos de área
  • Ingestión y manipulación de datos geoespaciales a gran escala para Deep Learning

5. Trabajo futuro

  • 🐛 Lista de tareas https://github.com/satellogic/telluric/issues
  • Mejor visualización (buen rendimiento para grandes volúmenes, más interactividad)
  • API aún más sencilla
  • Mejor manejo de casos límite (antimeridiano, polos)
  • Más operaciones geométricas avanzadas (buffering geodésico, centroides)
  • ...Sorpresas